Lecture 4
NHMRC Clinical Trials Centre, University of Sydney
gillian.heller@sydney.edu.au
Select the response distribution
Select covariates for all of the model parameters of the chosen distribution
Check model fit
Go back to 1.
y_i-\hat y_i
Suitably standardised, these are approx normally distributed if model assumptions are correct
This doesn’t extend to non-normal response distributions
We use quantile residuals (or probability integral transform (PIT) residuals)
F(y_i)\sim\mathcal{U}(0,1) (This is a well-known result.)
We fit a model and obtain parameter estimates (\hat\theta_{i1},\ldots,\hat\theta_{iK})
Compute
\hat u_i=F\left(y_i|\hat\theta_{i1},\ldots,\hat\theta_{iK}\right)
\hat r_i=\Phi^{-1}(\hat u_i)
The \hat r_i are the normalized quantile residuals, which we can
test for normality
use for identification of outlying observations
This means that randomized quantile residuals will be (slightly) different each time they are computed
set.seed())The worm plot is a detrended version of the QQ plot
Lack of fit in the tails is highlighted better
Can be split by covariate for identifying where lack of fit is occurring
Recall the full data set
We fit the following model
pb() = P-spline
The worm plot indicates a mostly well-fitting model.
The following divides age into 6 levels with approx equal numbers of observations:
number of missing points from plot= 1 out of 1222
number of missing points from plot= 3 out of 1209
number of missing points from plot= 1 out of 1217
number of missing points from plot= 4 out of 1215
number of missing points from plot= 0 out of 1215
number of missing points from plot= 1 out of 1216
The following specifies cutpoints for age:
Another model for BMI: Gamma (GA) distribution
k=2 i.e. AICFor the BCTo distribution, \mu is the median.
pred <- predict(m.BCTo, what="mu", type = "response")
dbbmi <- cbind(dbbmi, pred)
ggplot(dbbmi, aes(x=age, y=bmi)) +
geom_point(shape = 21, fill = "lightgray", color = "black", size = 0.5) +
geom_smooth(aes(x=age, y=pred), size=0.5, colour = "blue") +
theme_bw()